OpenCyto Analysis of HVTN080
The HVTN080 study flow cytometry data has been preprocessed (loaded in R from the raw FCS files and the FlowJo workspaces). The gating set is available here, and has been gated using OpenCyto.
We have gated the data using this OpenCyto template.
The manually gated data is available here, and contains the manual gates from the FlowJo workspaces.
The raw FCS files are available from flowrepository.org/id/FR-FCM-ZZ7U.
The data in the gating sets has already been compensated and transformed using the compensation and transformation information stored in the FlowJo workspaces.
The FlowJo workspaces for this data will be posted soon.
The code to perform the gating and extract the features is here.
Visualization of the gating hierarchy
We load the gated data and visualize the results.
suppressPackageStartupMessages(
{require(flowWorkspace)
require(flowViz)
require(knitr)
require(data.table)
require(MIMOSA)
require(GFMisc)
require(scales)
require(lme4)
require(contrast)
require(multcomp)
require(COMPASS)
})
opts_chunk$set(list(message=FALSE,warning=FALSE))
WORKDIR<-"/Users/gfinak/Dropbox/GoTeam/Projects/Paper-OpenCytoPipeline/HVTN080/reproducible/"
AUTO<-"./gs_auto_clean"
MANUAL<-"./gs_manual_clean"
setwd(WORKDIR)
auto<-load_gs(AUTO)
## loading R object...
## loading tree object...
## Done
manual<-load_gs(MANUAL)
## loading R object...
## loading tree object...
## Done
Here is a view of the gating hierarchy for the manually gated data.
plot(manual[[1]])
And here is the same for the OpenCyto gated data.
plot(auto[[1]])
Visualize gating layouts for manual and automated gating
We can view dotplots of each gate from selected samples for each of the OpenCyto and manually gated data.
First, the OpenCyto gated data:
plotGate(auto[[1]],scales=list(cex=c(1.2,1.2)),par.strip.text=list(cex=1.4),xlab=list(fontsize=12),ylab=list(fontsize=12),par.settings=list(gate.text=list(cex=1)),path=2)
Then the manually gated data:
plotGate(manual[[1]],scales=list(cex=c(1.2,1.2)),par.strip.text=list(cex=1.4),xlab=list(fontsize=12),ylab=list(fontsize=12),par.settings=list(gate.text=list(cex=1)),path=2)
Comparison of the perforin gate
We want to look at the variablity of the perforin channel in the OpenCyto analysed data.
#extract manually gated perforin for CD8+ for a single sample
#p1<-plotGate(manual["505325.fcs"],"8+/Perforin+",scales=list(cex=c(1.2,1.2)),par.strip.text=list(cex=1.1),xlab=list(fontsize=18),ylab=list(fontsize=18),par.settings=list(gate.text=list(cex=2)),isPath=FALSE,main="Manual Gate",arrange=FALSE)
#Extract the same gate from the OpenCyto gated data
set.seed(100)
s<-sample(1:length(auto),6)
p2<-plotGate(auto[s],"8+/Perforin+",scales=list(cex=c(1.1,1.1)),par.strip.text=list(cex=1.1),xlab=list(fontsize=12),ylab=list(fontsize=12),par.settings=list(gate.text=list(cex=1)),isPath=FALSE,default.y="<Alexa 680-A>",main="OpenCyto Gate",arrange=FALSE)
#generate a nice layout in a single panel
pdf(file="../../Submission/Figure4.pdf",width=8,height=5)
grid.arrange(p2,main="Perforin Expression")
dev.off()
## pdf
## 2
tiff(file="../../Submission/Figure4.tiff",width=8*150,height=5*150,res=150)
grid.arrange(p2,main="Perforin Expression")
dev.off()
## pdf
## 2
grid.arrange(p2,main="Perforin Expression")
Side-by-side comparison of OpenCyto and manual gating for other gates
We can generate any layout we want, even viewing plots for the same gate and sample side-by-side from the manual and openCyto analysis.
We can control the display projections of individual gates.
# Comparison of manual vs Autogating for ICS
proj=list("S"=c(x="FSC-H",y="FSC-A"),"L"=c(x="FSC-A",y="SSC-A"),"4+"=c(x="CD8",y="CD4"),"3+"=c(y="CD3",x="CD8"),"8+/57+"=c(x="CD57",y="Granzyme B"),"8+/Granzyme B+"=c(x="Perforin",y="Granzyme B"),"4+/IL2+"=c(x="IL2",y="IFNg"),"4+/IFNg+"=c(x="IL2",y="IFNg"),"8+/IL2+"=c(x="IL2",y="IFNg"),"8+/IFNg+"=c(x="IL2",y="IFNg"))
man<-plotGate(manual[["505325.fcs"]],c("S","Lv","L","3+","4+","8+","4+/IFNg+","4+/IL2+","8+/IFNg+","8+/IL2+","8+/Granzyme B+","8+/57+"),projections=proj,scales=list(cex=c(0.7,0.7)),par.strip.text=list(cex=0.7),xlab=list(fontsize=12),ylab=list(fontsize=12),par.settings=list(gate.text=list(cex=0)),path=1,arrange=FALSE,xbin=64,merge=TRUE)
aut<-plotGate(auto[["505325.fcs"]],c("S","Lv","L","3+","4+","8+","4+/IFNg+","4+/IL2+","8+/IFNg+","8+/IL2+","8+/Granzyme B+","8+/57+"),projections=proj,scales=list(cex=c(0.7,0.7)),par.strip.text=list(cex=0.7),xlab=list(fontsize=12),ylab=list(fontsize=12),par.settings=list(gate.text=list(cex=0)),path=1,arrange=FALSE,xbin=64,merge=TRUE)
gates<-list()
for(i in 1:9){
gates[[i]]<-c(man[[i]],aut[[i]])
}
grid.arrange(do.call(arrangeGrob,gates),ncol=1)
pdf(file="../../Submission/Figure2.pdf",width=8,height=6)
grid.arrange(do.call(arrangeGrob,gates),ncol=1)
dev.off()
## pdf
## 2
tiff(file="../../Submission/Figure2.tiff",width=8*150,height=6*150,res=150)
grid.arrange(do.call(arrangeGrob,gates),ncol=1)
dev.off()
## pdf
## 2
Evaluation of vaccine-specific effects in Perforin-positive subsets of cells
Since OpenCyto was able to accurately gate perforin, we’d like to see if we can detect any vaccine-specific perforin cell subsets.
We extract cell counts from the CD4 and CD8 T-cell subsets that express perforin in combination with other cytokines, using COMPASS.
First, we extract the counts from the automated gates using COMPASS. We focus on cells expressing any cytokine together with perforin, and test whether the proportion of pre vs post-vaccine cells differs by treatment group (Wilcoxon test) within the CD4 and CD8 subsets.
#extract the population statistics
auto_stats<-getPopStats(auto,statistic="count")
if(!file.exists("./CD8_CompassContainer.rds")){
cd8_auto<-COMPASSContainerFromGatingSet(gs=auto,node="8\\+$",individual_id="PTID",sample_id="name",markers=c("IFNg","TNFa","IL2","Granzyme B","Perforin"))
saveRDS(cd8_auto,file="/Users/gfinak/Dropbox/GoTeam/Projects/Paper-OpenCytoPipeline/HVTN080/CD8_CompassContainer.rds")
}else{
cd8_auto<-readRDS(file="./CD8_CompassContainer.rds")
}
if(!file.exists("./CD4_CompassContainer.rds")){
cd4_auto<-COMPASSContainerFromGatingSet(gs=auto[pData(auto)$Stim%in%c("POL-1-PTEG","negctrl 1","negctrl2 ")],node="4\\+$",individual_id="PTID",sample_id="name",markers=c("IFNg","TNFa","IL2","Granzyme B","Perforin"))
saveRDS(cd4_auto,file="/Users/gfinak/Dropbox/GoTeam/Projects/Paper-OpenCytoPipeline/HVTN080/CD4_CompassContainer.rds")
}else{
cd4_auto<-readRDS("./CD4_CompassContainer.rds")
}
#Plots the proprotions for a given combination for CD8
.plotProps<-function(cd8,title,combination){
#Compute cell counts for the given combination
C<-COMPASS:::CellCounts(cd8,combinations=combination)
tofit<-cbind(cd8$meta,cbind(cd8$counts-C,C))
colnames(tofit)[c(25,26)]<-c("NSUB","CYTNUM")
tofit$Stim<-factor(tofit$Stim)
#Aggregate negative controls
levels(tofit$Stim)<-c("ENV-1-PTEG","GAG-1-PTEG","negctrl","negctrl","POL-1-PTEG")
#relevel the visit variable
levels(tofit$VISITNO)<-c("Pre-vaccine","Post-vaccine")
tofit$rx_code<-factor(tofit$rx_code)
#relevel the treatment variable
levels(tofit$rx_code)<-c("Placebo","Pennvax B", "Pennvax B + IL12 DNA")
tofit2<-ddply(tofit,.(Stim,VISITNO,PTID),function(x)with(x,data.frame(rx_code=unique(rx_code),Stim=unique(Stim),CYTNUM=sum(CYTNUM),NSUB=sum(NSUB))))
p1<-ggplot(subset(tofit2,Stim%in%"POL-1-PTEG"))+geom_boxplot(aes(x=rx_code,y=CYTNUM/(NSUB+CYTNUM)))+facet_grid(Stim~VISITNO)+theme_bw()+scale_y_continuous("Empirical Proportion of Vaccine\nResponsive Cells")+scale_x_discrete("Treatment Group")+coord_trans(ytrans=log1p_trans())
E<-ConstructMIMOSAExpressionSet(tofit2,reference=Stim=="negctrl",measure.columns=c("NSUB","CYTNUM"),other.annotations=c("rx_code","VISITNO","Stim","PTID"),default.cast.formula=component~PTID+Stim+VISITNO+rx_code,.variables=.(PTID,rx_code,VISITNO))
fitted<-MIMOSA(NSUB+CYTNUM~rx_code+PTID+RefTreat+Stim+VISITNO,E,ref=RefTreat%in%"Reference"&Stim%in%"POL-1-PTEG",subset=RefTreat%in%"Treatment"&Stim%in%"POL-1-PTEG",method="mcmc",pXi=1,EXPRATE=1e-6,tune=500)
#linear model contrast pr vs post vaccine within treatment groups
lmdata<-cbind(with(data.frame(countsTable(fitted,proportion=TRUE)),data.frame(proportion=CYTNUM-CYTNUM_REF)),"Pr.Response"=getZ(fitted)[,2],pData(fitted))
KK<-GFMisc:::optimAsinhCofactor(subset(lmdata,rx_code%in%"Placebo")$proportion)
lmdata$prop_trans<-asinh(lmdata$proportion*KK)
lmfit<-lm(prop_trans~rx_code*VISITNO,data=lmdata)
contr<-contrast(lmfit
,a=list(rx_code=c("Pennvax B","Pennvax B + IL12 DNA"),VISITNO="Post-vaccine"),
b=list(rx_code=c("Pennvax B","Pennvax B + IL12 DNA"),VISITNO="Pre-vaccine"))
fit1<-(lmer(prop_trans~rx_code*VISITNO+(1|PTID),data=lmdata))
sumry<-data.frame(p.value=summary(glht(fit1,linfct=contr$X,alternative="greater"))$test$pvalues,rx_code=c("Pennvax B","Pennvax B + IL12 DNA"))
lmdata<-data.table(lmdata)
setkey(lmdata,PTID,VISITNO)
#extract the data and plot
lmdata_plot<-lmdata[,list(difference=proportion[2]-proportion[1]),list(rx_code,PTID)]
p<-ggplot(subset(lmdata_plot,!rx_code%in%"Placebo"),aes(x=rx_code))+geom_boxplot(aes(y=difference))+theme_bw()+geom_text(aes(y=3e-4,label=paste0("p = ",signif(p.value,3))),data=sumry)+scale_y_continuous("Post-Vaccine - Pre-Vaccine")+scale_x_discrete("Treatment Group")
return(p)
}
#Plots the proportion for a given combination for CD4
.plotPropsCD4<-function(cd8,title,combination){
C<-COMPASS:::CellCounts(cd8,combinations=combination)
tofit<-cbind(cd8$meta,cbind(cd8$counts-C,C))
colnames(tofit)[c(25,26)]<-c("NSUB","CYTNUM")
tofit$Stim<-factor(tofit$Stim)
#Aggregate negative controls
levels(tofit$Stim)<-c("negctrl","POL-1-PTEG")
#relevel the pre-post vaccine variable
levels(tofit$VISITNO)<-c("Pre-vaccine","Post-vaccine")
tofit$rx_code<-factor(tofit$rx_code)
#relevel the treatment groups
levels(tofit$rx_code)<-c("Placebo","Pennvax B", "Pennvax B + IL12 DNA")
#some data reshaping
tofit2<-ddply(tofit,.(Stim,VISITNO,PTID),function(x)with(x,data.frame(rx_code=unique(rx_code),Stim=unique(Stim),CYTNUM=sum(CYTNUM),NSUB=sum(NSUB))))
E<-ConstructMIMOSAExpressionSet(tofit2,reference=Stim=="negctrl",measure.columns=c("NSUB","CYTNUM"),other.annotations=c("rx_code","VISITNO","Stim","PTID"),default.cast.formula=component~PTID+Stim+VISITNO+rx_code,.variables=.(PTID,rx_code,VISITNO))
# We don't need to fit a MIMOSA model here, I just want to extract the cell counts with the metadata.
fitted<-MIMOSA(NSUB+CYTNUM~rx_code+PTID+RefTreat+Stim+VISITNO,E,ref=RefTreat%in%"Reference"&Stim%in%"POL-1-PTEG",subset=RefTreat%in%"Treatment"&Stim%in%"POL-1-PTEG",method="mcmc",pXi=1,EXPRATE=1e-6,tune=500)
#linear model testing pre vs post-vaccine within each treatment group
lmdata<-cbind(with(data.frame(countsTable(fitted,proportion=TRUE)),data.frame(proportion=CYTNUM-CYTNUM_REF)),"Pr.Response"=getZ(fitted)[,2],pData(fitted))
KK<-GFMisc:::optimAsinhCofactor(subset(lmdata,rx_code%in%"Placebo")$proportion)
lmdata$prop_trans<-asinh(lmdata$proportion*KK)
lmfit<-lm(prop_trans~rx_code*VISITNO,data=lmdata)
contr<-contrast(lmfit
,a=list(rx_code=c("Pennvax B","Pennvax B + IL12 DNA"),VISITNO="Post-vaccine"),
b=list(rx_code=c("Pennvax B","Pennvax B + IL12 DNA"),VISITNO="Pre-vaccine"))
fit1<-(lmer(prop_trans~rx_code*VISITNO+(1|PTID),data=lmdata))
sumry<-data.frame(p.value=summary(glht(fit1,linfct=contr$X,alternative="greater"))$test$pvalues,rx_code=c("Pennvax B","Pennvax B + IL12 DNA"))
lmdata<-data.table(lmdata)
setkey(lmdata,PTID,VISITNO)
lmdata_plot<-lmdata[,list(difference=proportion[2]-proportion[1]),list(rx_code,PTID)]
#Pull out the results and plot
p<-ggplot(subset(lmdata_plot,!rx_code%in%"Placebo"),aes(x=rx_code))+geom_boxplot(aes(y=difference))+theme_bw()+geom_text(aes(y=3e-4,label=paste0("p = ",signif(p.value,3))),data=sumry)+scale_y_continuous("Post-Vaccine - Pre-Vaccine")+scale_x_discrete("Treatment Group")
return(p)
}
#proportion of cells expressing any cytokine and perforin, by treatment group
cap<-capture.output(anycyt<-.plotProps(cd8=cd8_auto,title="POL-1-PTEG Stimulation of\nCD8+ T-cells for (IL2 or IFNg or TNFa) AND Perforin",combination="(IL2|IFNg|TNFa)&Perforin"))
cap<-capture.output(cd4_anycyt<-.plotPropsCD4(cd8=cd4_auto,title="POL-1-PTEG Stimulation of CD4+ T-cells\nfor (IL2 or IFNg) AND Perforin",combination="(IL2|IFNg)&Perforin"))
grid.arrange(arrangeGrob(arrangeGrob(anycyt,main="CD8 (Any Cytokine & Perforin)"),arrangeGrob(cd4_anycyt,main="CD4 (Any Cytokine & Perforin)"),nrow=1))
SessionInfo for reproducibility
sessionInfo()
## R Under development (unstable) (2014-02-06 r64933)
## Platform: x86_64-apple-darwin13.0.0 (64-bit)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] splines parallel grid stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] sandwich_2.3-0 COMPASS_1.1.5-1 multcomp_1.3-2
## [4] TH.data_1.0-3 mvtnorm_0.9-9997 contrast_0.19
## [7] rms_4.1-3 SparseM_1.03 Hmisc_3.14-0
## [10] Formula_1.1-1 survival_2.37-7 lme4_1.0-6
## [13] Matrix_1.1-2 scales_0.2.3 GFMisc_1.0
## [16] MIMOSA_0.99.5 ggplot2_0.9.3.1 Biobase_2.23.4
## [19] BiocGenerics_0.9.3 reshape_0.8.4 plyr_1.8
## [22] MASS_7.3-29 data.table_1.8.11 knitr_1.5.22
## [25] flowWorkspace_3.9.65 gridExtra_0.9.1 ncdfFlow_2.9.19
## [28] flowViz_1.27.16 lattice_0.20-27 flowCore_1.29.23
##
## loaded via a namespace (and not attached):
## [1] abind_1.4-0 clue_0.3-47 cluster_1.14.4
## [4] coda_0.16-1 colorspace_1.2-4 corpcor_1.6.6
## [7] DEoptimR_1.0-1 dichromat_2.0-0 digest_0.6.4
## [10] evaluate_0.5.1 feature_1.2.10 formatR_0.10
## [13] graph_1.41.2 gtable_0.1.2 hexbin_1.26.3
## [16] IDPmisc_1.1.17 KernSmooth_2.23-10 ks_1.8.13
## [19] labeling_0.2 latticeExtra_0.6-26 MCMCpack_1.3-3
## [22] minqa_1.2.3 modeest_2.1 munsell_0.4.2
## [25] nlme_3.1-113 pcaPP_1.9-49 pracma_1.6.4
## [28] proto_0.3-10 RColorBrewer_1.0-5 Rcpp_0.11.0.2
## [31] reshape2_1.3.0.99 Rgraphviz_2.7.0 rmarkdown_0.1.5
## [34] robustbase_0.90-2 rrcov_1.3-4 RSVGTipsDevice_1.0-4
## [37] stats4_3.1.0 stringr_0.6.2 testthat_0.8.1
## [40] tools_3.1.0 XML_3.98-1.1 yaml_2.1.10
## [43] zlibbioc_1.9.0 zoo_1.7-11